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Abstract 

We show theoretically that the hypothesis of criticality as a theory 
of long-range fluctuation in the human brain may be distinguished from 
competing explanations on the basis of macroscopic neuronal signals such 
as the EEG, using novel theory of narrowband amplitude time-series at 
criticality. Our theory predicts the division of critical activity into meta¬ 
universality classes. As a consequence we provide strong evidence for 
criticality in the human brain on the basis of the EEG. 

Testing for the presence of critical power-law avalanche dynamics (CPLAD) 
is pivotal in establishing that a dynamical system is poised close to a critical 
point. Such an analysis assumes that periods of activity are separated from 
periods of inactivity. For the intact brain of awake human subjects, which 
has been hypothesized to be critical [T|, the assumption of separation between 
activity and inactivity is never met, since the system is continuously active. This 
lack of separation makes unclear whether the results of existing experiments 
confirming criticality in animals in vitro Hum and in vivo HE] generalize to 
the intact brain of awake human subjects; in particular it has been argued that 
the scale-free form of the power-spectra of LFP, EEG and MEG data is due 
not to CPLAD but to passive filtering (PF) through the extracellular media 
While several papers have claimed to test the hypothesis of criticality 
from large scale-brain signals of awake human subjects [5] the authors only 
find CPLAD over short time-scales (< 1 second) since the continuous nature 
of the data prevented testing of criticality as a theory of neuronal fluctuations 
over longer time-scales. Quantifying potential CPLAD from such non-invasive 
human recordings is challenging since a superposition of overlapping avalanches 
may lead to incorrect estimation of size and duration (critical) exponents (this 
is related to the under sampling problem 0) or indeed of whether CPLAD are 
present at all. Interestingly, in this letter we show that in order to quantify and 
distinguish CPLAD from PF, it is not necessary to define individual avalanches; 
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instead we show that is is sufficient to analyse the properties of continuous 
neuronal recordings. This approach allows us to test the hypothesis of criticality 
as a theory of neuronal variability over large time-scales. 

We use the facts that periods of activity of local neural networks are sepa¬ 
rated by periods of inactivity, and macroscopic brain signals, such as the LFP, 
EEG and MEG, measure the linear superposition X'(t) of activity of numer¬ 
ous local neural networks f5J; this means that the macroscopic brain signal is 
continuous, although the activity of local networks is not. Thus: 


T q. 

X'{t) = Y Y. h a<i as,i 

s =1 i =1 



(i) 


a s,i \tY) denotes activity at time t of the * th local neural network L s i which 
begins at time s, after a period of inactivity and lasts for L Sti time steps. We 
adopt the conventions that a Sj ,;(f) is normalized to unit standard deviation and 
is 0 for t £ (—oo,0) U (l,oo). Thus, h s ^ denotes the average height of the 
time-course of activity of the i th local neural network which begins at time s. 

PF and CPLAD differ on how h s ^ and L St i are distributed. According to the 
PF hypothesis their distributions decay faster than power-laws. For simplicity, 
therefore, we summarize the PF as claiming exponential distributions: p(L) ~ 
e~ AL and p{h) ~ e~ Bh . Moreover the PF hypothesis states power-spectra 
of macroscopic brain signals X(t) are scale-free because they reflect a filtered 
version of X'(t): 


X(t) = yF(u)X'(t-u) (2) 

u =0 

F(u) is a linear filter due to the extracellular media which yields a signal with 
power-law power spectrum from white noise input El El- 

On the other hand, the CPLAD hypothesis states that macroscopic signals 
X(t) reflect X'(t) directly so that X(t) = X'(t) and that we have power law 
distributions for h s ^ and L s i : L ~ L~ a and h ~ ZA These power-law dis¬ 
tributions explain the power-law form of the power-spectrum m- In addition 
the CPLAD hypothesis states that each a a ^(t) is an independent and identical 
sample from a single stochastic process, which we call a(t) [Tl] . 

We now present theory which makes predictions for the CPLAD hypothesis 
which provably do not hold for the PF hypothesis. As well as considering the 
raw signal X(t) we also consider the amplitude of a narrowband filtered version 
of X(t), which we denote g u (X(t)). Thus let /„(•) be a linear narrowband filter 
with pass band w G [w — Au],w + Aw], and H(-) the Hilbert transform then: 

gUX(t)) = | MX(t)) + iH(f u (X(t))) | 

The theory depends on two measures of long-term variation evaluated on 
X(t) and g 0J (X(t)): the Hurst exponent and the detrended cross correlation 
coefficient. The Hurst exponent H, of a stationary process Y ( t ) may be defined 
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by the scaling of the auto covariance function. For H G (0,1), Y(t) is said to 
be long-range temporally correlated (LRTC): 

E(Y(i + s)Y(t)) - E (Y{t)) 2 ~ s 2H ~ 2 


For auto covariances decaying faster than s _1 , one defines H = 1/2 and Y(t) 
is not LRTC. From now on we distinguish the Hurst exponents of X(t) and 
g bJ (X(t)) as H raw and H“ mp . Typically in neuroscientific applications Hurst 
exponents are measured with Detrended Fluctuation Analysis (DFA) [12] , 

The detrended cross correlation coefficient [13] pDCCA{n , Y\^Y 2 ) is a measure 
of correlation between two time-series Y\{t) and Y 2 {t) at a time-scale n, which 
is invariant to non-stationary trends of a fixed polynomial degree. Let F Y Yl ( n ) 
be the detrended variance of Y\ where each window of size n is detrended, 
as computed for DFA and F Y Y2 (n), by analogy, the detrended covariance, as 
computed for DCCA El; then in analogy to the Pearson correlation coefficient, 
one defines: 


PDCCA(n , Yi,Y 2 ) 



(n) 


Our first result is that for the X(t) of the PF hypothesis, Equation |5]), H£ mp = 
0.5 and, when uq ^ w 2 , PDCCA{n, g UI1 (X(t)), g(uj 2 ){X(t))) ->0asrt->oo. This 
can be seen by splitting Equation @ into activity which lasts longer than some 
value L' and activity which has duration shorter than or equal to L 1 : 


X'{t) 




(3) 


Since the distribution of L decays exponentially, the vari- ance of the left hand 
term dominates. Therefore: 


L s i<L' ' \ S A / / 


Since time-points in this expression spaced more than L' points apart are inde¬ 
pendent we also have that time-points of g u {X'(t)) spaced more than L' points 
apart are independent, so that H^ mp = 1/2 and PDCCA(n,u}i,uj 2 ) —> oo as 
n —> oo (same asymptotic properties as white noise); the reweighting of fre¬ 
quency induced by the passive filtering does not change the narrowband prop¬ 
erties. See Proposition 2.1 of the supplement. 

We now consider the CPLAD hypothesis. The behaviours we derive for 
Hamp an d PDCCA{n, (X(t)), g U2 (X(t))) depend on the exponents a and /?. 
We find that there are four regions of parameters with qualitatively differing 
behaviours which we term meta-universality classes. We assume that all power- 
law distributions are cut off at a lifetime L c which is proportional to the size of 
the system and, for simplicity, that at each time-point a fixed number q s = q of 
avalanches begin. 
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(MU1) a < 2 


For a < 2 the number of avalanches active at time t is: 


^ ^ /({Lt — s,i\I J t—s,i ^ ^ 


l 


L, 


u 


L. 


L~ a dL 



This implies the number of avalanches active at any given time is unbounded in 
the system size. Applying the Central Limit Theorem, this implies that X(t) is 
a Gaussian process with power-law autocorrelation (see supplement for details). 
Since Gaussian processes are uniquely defined by their second order properties 
and may be generated by filtering white noise he then in analogy to the results 
for the PF hypothesis H£ mp = 1/2 and pDCCA{n, uj±, 0 J 2 ) —» 0 as n —> 00 . See 
Proposition 2.4. 

(MU2) a > 2 and a < (3 + 3 

For a > 2, the probability that an avalanche is active with duration greater 
than L' is: 



Ls,!>L' 


Thus the probability that large avalanches occur simultaneously is negligible. 
Moreover we have, for an exponent (3', a scaling relation for large L: 


fui{L^a{t/L)) ~ if f u (a){t/L) 


Here f w (a)(t/L ) is understood as position t/L of a(t) filtered in the narrowband 
around w. We are able to derive this exponent theoretically (Proposition 2.7) 
and find that f3' = (3/2. This implies that for a < (3 + 3 (but not otherwise) the 
variance of the large narrowband filtered avalanches dominates so that: 



Since the avalanches on the right hand side of this relation do not overlap then 
we have that: 



Given this representation of the amplitude as a simple sum of amplitudes of in¬ 
dividual avalanches, standard techniques may then be applied to approximating 
the Hurst exponent H% m, and we find: 




Similarly we find: 


H? aw =P~ a/2 + 2 

> 1/2 


Moreover, assuming that the integrals J/' g u {a{t))dt exist, then the separation of 
large avalanches make it simple to derive that pDCCA(n, uq, uq) —t 1 as n — > 00 . 
See Propositions 2.8 and 2.9. 

(MU3) a > 2 and a < 2/3 + 3 

For frequencies with 1 /uj /$> L c and L < L c : 

f u (L 0 )~LPf u (a){t/L) 

This is because relative to the time-scale of the filter, a(t/L) may be treated 
as a delta function, which weights all frequencies equally. Therefore applying a 
similar argument as for MU2 we find that for uq , uq — > 0, pdcca {? i , uq, uq) —► 
lasn—> 00 . Moreover, applying the results for MU2 we have: 

Km P = 1/2 
H? aw =/3^a/2 + 2 
> 1/2 


See Proposition 2.10. 

(MU4) a > 2/3 + 3 and a > 2 

Since these universality classes have the shortest tails in their critical distri¬ 
butions, we may apply identical methods as applied to the PF hypothesis which 
show that H£ mp = 1/2 and pdcca{ti, uq, UI 2 ) —> 0 as n —> 00 when uq ^ oq. 
See Proposition 2.11. 

The results of the theory for the CPLAD model are summarized in Figure [T] 

We now present tests of these predictions in simulations. We first tested 
the predictions for the PF hypothesis by generating a LRTC process by linear 
filtering of white noise. The results are displayed in the bottom panel of Figure [2] 
and Figure 2 of the supplement. The results confirm that H^ mp = 1/2 and 
PDCC J 4(n,wi,w 2 ) -t 0 . 

We then tested the predictions for the CPLAD hypothesis, modeling the 
activity of local networks by: 


a(t) = b(t) + c(t)e(t ) 

b(t) is the average avalanche shape, c(t) the shape in the variance profile and e(f) 
a colored noise with the spectrum known for a critical system V[uj\ ~ p~6] 

(see supplement for details). Two examples in MU2 and MU1 are displayed 
in the top two rows of Figure [2] 

The results of a simulation for all critical parameters are displayed in Fig¬ 
ure [ 3 ] We find good agreement be- tween the meta-universality classes derived 
and the simulation results. (See supplement sec. 3.3 for details.) 
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Behaviour over critical parameters 



Figure 1: Division of critical exponents a and /3 into rneta- universality classes. 
The figure displays the range of qualitative behaviours we predict with our 
theory. Areas marked in green display no LRTC behaviour in sub-bands or 
DCCA correlations between sub-bands. Areas in red display LRTC and/or cross 
correlations between amplitudes of sub-bands (H% mp = 1/2, Pdcca{p) = 0 for 
large n ). 


6 







1 


0.5 






n - time scale 



Passive Filtering Model 



Figure 2: Right: Sample paths from the PF and CPLAD models. In each of the 
three cases the x-axis denotes time and the y-axis number of act i vat ions. In each 
case the middle trace denotes X(t) and the top and bottom g u . ( X(t )). Left: 
DCCA correlation coefficients pdcca(ji ) and Hurst exponents H^ mp . 



Figure 3: DCCA correlation coefficients and Hurst exponent for the simulated 
CPLAD model. The x-axis denotes the exponent a and the y-axis denotes (3. 
The black lines denotes the transitions in meta-universality class according to 
the theory. 
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Figure 4: Results of data analysis of human EEG. The frequency ranges analysed 
i = 1, 2, 3 are 35-40Hz, 60-65Hz and 72-77 Hz respectively, which are displayed 
as superscripts in the plots. 


Finally we tested the PF and CPLAD hypotheses by estimating H^ mp , H raw 
and pDCCA(n,uii,ui 2 ) on EEG time-series (see supplement for preprocessing). 
Important is that we analyse 3 frequency ranges without oscillations (no local 
maximum in power-spectrum); the aim was to restrict analysis to activity cor¬ 
responding to the 1 // 7 shape of the power-spectra. Given that the data were 
sampled at 200Hz, and that lower frequencies require far larger window sizes for 
analysis, we chose 3 frequencies above the beta range, taking care to exclude 
the 50Hz line noise. 

The results of our analysis over all 7 subjects are displayed in Figure |4j 
For 244 of 261 signals, the DFA estimate of H£ mp was higher than 1/2, clearly 
indicating the presence of LRTC (median = 0.61, 5th and 95th per¬ 

centiles 0.49 and 0.85, p j 0.0001). Likewise, 238 of 261 DCCA correlation 
coefficients pDCCA(n,uJi,u) 2 ) at the highest scale n between frequency ranges 
were positive (jp -C 0.0001). We found moreover that the Pdcca values at 
the highest scale and H^ mp measured from the same neural data were highly 
correlated (p <C 0.0001) and H£ mp values in distinct frequency ranges were 
highly correlated (p <C 0.0001) (likewise for pjjCCA(n)). Finally we found that 
the pDCCA{n ) and H^ mp values were not significantly correlated with H rau , 
(p > 0.05). These results confirm the CPLAD hypothesis, in particular in 
agreement with predictions for MU2. 

Thus these findings provide strong evidence for criticality in the human brain 
by confirming the hypothesis of CPLAD. These results cannot be explained by 
the PF hypothesis; passive filtering of neural signals may generate a power- 
law spectrum but will not induce LRTC of narrowband amplitudes or DCCA 
correlations between narrowband amplitudes. Moreover, we note that the failure 
to detect avalanches in the experiments of HE] may be explained by criticality 
in the first meta-universality class where we have shown that even the largest 
avalanches are not discernible. Thus it is possible in these experiments that 
the systems under study were critical despite failure to demonstrate CPLAD 
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A Estimators 

A.l Detrended Fluctuation Analysis 

Detrended Fluctuation Analysis (DFA) [T2] is a methodology for the estimation 
of the Hurst exponent of a (possibly non-stationary) time-series. Its advantage 
over covariance analysis or analysis of the power-spectrum are its robustness 
to trends contaminating the empirical time-series and its desirable convergence 
properties El- 

The steps involved in DFA are as follows. First one forms the aggregate sum 
of the empirical time-series X(t): 


t 

x (t) = J^ x (i) 

i—1 

(From now on whenever we refer to x(t) we mean the time-series obtained 
from X(t ) by way of this operation.) Analysis of the fluctuations in X{t) may 
then be performed by measuring the variance of x{t ) in windows of varying size n 

after detrending, i.e., x{t) is split into windows of length n, ri 1 ',..., Xn \ ..., ir^ A 
and the average variance after detrending the data of x{t) in these windows is 
formed; i.e. let Pd be the operator which generates the mean-squares estimate of 
the polynomial fit of degree d, then the DFA coefficients or detrended variances 
of degree d are: 
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Crucially, in the limit of data the slope of log(Fp FA (n)) against log(n) converges 
to H. Thus X(t) is power-law correlated LRTC if and only if the estimate of H , 
H , converges to a number greater than 0.5 in the limit of data. We note here 
that there are numerous methods for the Hurst exponent; these include wavelet 
estimators ilBj, log-periodogram based methods [15], among others [20]. We 
chose DFA since it is standard in the physics and neuroimaging literature, and 
yields competitive estimates [211 . 

A.2 Detrended Cross-Correlation Analysis 

In precise analogy to DFA, Podobnik and Stanley propose Detrended Cross- 
Correlation Analysis (DCCA) [2], an extension of DFA to two time-series, by 
considering the detrended covariance: 

FlccAin) = j . n (( Xl )™ ) ~ p d (OOn 5 )) ((*2 )n ] - Pd ((*2)£°)) 

Thus DCCA generalizes DFA in the sense that if X± = X 2 then F^ CCA (n) = 
P DFA(. n )- To simplify the fact that we will need to consider the DFA coefficients 
of A'i and X 2 simultaneously, we will also refer to the DCCA coefficients as: 

Fx 1 ,x 2 ( n ) = P D CCA ( n ) 

Thus, DCCA quantifies the behaviour of the covariance between A'i and X 2 
over a range of time-scales given by n. In analogy to the Pearson correlation 
coefficient, we may consider the detrended cross correlation coeffcient H3: 


pDGCA{n , X 1 ,X 2 ) 



(n) 


Pdgca quantifies the correlation between X± and X 2 over a range of time- 
scales. Note that while the machinery involved in the estimation of pdcca are 
more complex than for Pearson correlation, both coefficients estimate the same 
quantity for stationary time-series, not contaminated by trends. Thus Pdcca 
generalizes the Pearson correlation coefficient. Applicability to non-stationary 
time-series is particularly important for the neural data analysis. Whenever the 
context allows for no ambiguity we abbreviate pdcca(p , 1 A'i, X 2 ) to Pdcca{ti)- 


B Theory of the Avalanche Process 

We derive here properties which hold for the CPLAD model but not for the 
PF model. The theory focuses on the Hurst exponent of A (t) which we denote 
H raw , the Hurst exponent of g u (X(t)), which we denote H^ mp and the DCCA 
correlation coefficient pDCCA{n, g UJl (X(t)), g U2 (X(t))). Throughout we assume 
that all power-law distributions are cutoff at a lifetime L c , that q s = q is fixed 
and deterministic for all s and that t is greater than L c to ensure stationarity. 
First we consider the PF model. 
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Proposition B.l. For the PF model ¥.(p DC CA(n,g Ul (X(t)),g U2 {X(t)))) -> 0 
for n —)■ oo and H^ mp = 1/2. 

Proof. For the passive filtering model, the right hand term of Equation (3) of 
the main paper dominates. That is: 

Z(t) = 'y ' h S ’iQ, s ,i 

L s ,i<L' 



Time points of g u (Z(t)) greater than L' apart are independent so H(f mp = 1/2. 
Moreover E (pdcca{u, (.X(t )), g U2 (X(t)))) — > 0 since the numerator of the 
DCCA correlation coefficient is given by the DCCA coefficients F'f^ x 2 ( n )- 
These behave like the coefficients of white noise for large window sizes and 
therefore have expectation 0 in the limit. D 

We now turn to the CPLAD model, to which all of the following propositions 
apply. We first aim at understanding the scaling of the number of active large 
avalanches: 

Proposition B.2. The number of active avalanches is asymptotically ~ L^.~ a 
for large L c and a < 2 and of constant order for a > 2. 

Proof. The number of active avalanches at time t is equal to: 


J-'c 

4f{Lt-t',i\Lt~t',i > f} 

t '=0 



~ (1 - a)(ic) 1-Q - (1 - of) / {t'f- a dt' 

Jo 

= (1 - a)(L c y~ a - (1 - a)(2 - a)(L 2 ~ a - 1 ) 

j ~ Lc~ a for a < 2 
1 = 0(1) for a > 2 


□ 


Proposition B.3. X(t) is approximately Gaussian for large L c and for a < 2. 


Proof. At any time point we have a number k ~ L%T a (Proposition B.2) of 
avalanches which commenced each at times Sj = 1,..., k with lengths Lj. Then 
we may write: 


*(*) = !>< 


t - Si 
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For k large, which may be assumed by Proposition B.2 after normalization and 
centering we apply the Lyapunov central limit condition [22] to: 


X(t) = 


1 


VkL 0 c +{ 


a+l)/2 




3 =1 


t - Si 


where p is the mean value of aft) over the interval [0,1]. Let s\ = kL^ +a+1 , 
then for the Lyapunov condition to apply we need to show, for some 6 > 0, that: 


Thus: 


i k 

-2+1 E 


E 


3= 1 


U 



2+5 

-+ 0 


2+5 


E E 


3=1 






2+5 


1 

1+5/2 ^(/3+ a / 2 +l/2)(2+<5) 

1 

f : l+5/2]r / ^+a/2+l/2)(2+5) 

L a c +1 

k 5/2 L ( c a / 2 + 1 / 2 )( 2 + S ) 


Z L f+ s)p 

2=1 

fc ^/3(2+5)+a+l 


< k s t 2 L s J 2 


This completes the proof for univariate Gaussianity. The proof that con¬ 
vergence is to a multivariate Gaussian process is a simple extension via the 
Cramer-Wold device [22] . 

□ 


Proposition B.4. For a < 2, E(p DC cA(n,g Ul (X(t)),g U2 (X(t)))) -+ 0 for 
n -+ oo and H“ mp = 1/2. 


Proof. By Proposition B.3 X(t ) is a Gaussian process, which is unicpiely defined 
by its first two moments. Over long time-scales, therefore, X(t) behaves like 
fractional Brownian noise [25] . One approach to generating fractional Brownian 
noise is by long-range filtering of a Gaussian white noise m- This brings us 
exactly into the situation of Proposition |B.1| □ 


Proposition B.5. The probability that avalanches of length L' occur simulta¬ 
neously vanishes for large L' when a > 2. 


Proof. In analogy to Proposition B.2 equal to L' is the time spent in large 
avalanches divided by the total recorded time: 
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<3®(1(z/,oo)(£)) 


1 

T 


E L < 

Li>L' 



L 1 ~ a dL 


L ,2 ~ a 


Thus the probability that two avalanches of length greater than L' occur 
->• 0 . □ 

Proposition B.6. The Hurst exponent H raw ofX(t) satisfies the approximate 
relation: 


/3 - § + 2 if t < 3 
1/2 ifr> 3 


Proof. Following m we approximate avalanches by a box function: 


a(t) 


1 for t £ [ 0 , 1 ] 
0 otherwise 


(5) 


The authors show that in this case the autocorrelation function r[t) = 
E(X(s)A'(s + t)) — E(AT(s )) 2 satisfies: 

/»oo />oo 

r(t) oc / (L-\t\) L~ a L 2l3 dL 

J\t\ Jo 

~ \t\ 2/3 ~ a+2 

Using the fact that if the autocorrelation function scales as r(t) ~ t -7 , then the 
Hurst exponent and 7 are related as 7 = 2 — 2 H [24], so that: 

Ct 

H r aw = ft ^ + 2 ( 6 ) 

□ 

We adopt the convention that by g u (a)(t) we refer to the t th time-point of 
g^iait)). At criticality we have the relationship: 

gui{L?a(t/Li)) ~ Lf g UJ (a)(t/L i ) 


This scaling relationship is illustrated in Figure[5] The left hand panel illustrates 
avalanches whose heights scale according to Lr, with /3 = 1. The right hand 
panel displays g UJ (a(t/L)). One sees that the heights of these filtered avalanches 
scale more slowly than the original avalanches. 
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Raw Avalanches 


Filtered avalanches 




Figure 5: The figure illustrates the difference in scaling between the raw 
avalanches and their filtered amplitudes, when f3 = 1. The left hand panel dis¬ 
plays avalanches L^a{t/L) with log-spaced lifetimes L between 400 and 7000. 
The right hand panel displays the narrowband amplitudes of these avalanches 
L 13 g ul (a){t/L). Since (3 > f3' we see that the narrowband amplitudes scale less 
steeply than the raw avalanches. See Section |B] 


We require an additional known relation between critical exponents. Let S 
be the size of an avalanche, i.e. the total activity occurring over the course of 
the avalanche: S = L^a(t/L)dt. Then, at criticality: 

S~ S~ T 

Then we have the following relation [25] between critical exponents: 


a — 1 

T — 1 


= /3 + l 


If we define a' = a then we may also define t' by Equation Q so that: 

m-r+i 


(7) 

( 8 ) 


The size distribution is important for our subsequent analyses because whether 
the asymptotic properties of the process X (t) are dominated by the large avalanches 
depends on whether the size distributions have infinite variance or not. 
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Proposition B.7. f3' = /3/2. 

Proof. We have that V[uj\ ~ This implies that the standard deviation 

of /i w ,L 4 u (a(<)) scales according to L -/3 / 2 . This is because: 


And therefore: 


var(/ w 

, At 0 


pLJ-\-A(jJ 
/ w-Aw 


■p[a/]<iu/ 


1 

oT+i 


vai'(f Luj ,LAuj (a(t))) 


p Luj L Alo 

' Luj—LAuj 
p Luj-\~L Acu 


V[uj']duj' 


1 


' Luj—LAlo ^ 
puj-\-Aoj 


'/3+1 

1 


o—Auj mr +i 


duj' 

Ldui' 


■j-pVax(f„ A „(a(t))) 


Therefore: 

Std(fLu:,LAcj(a(t)))) ~ A _,3/2 

and: 

fu.Au,( L i ai(t/Li )) ~ Lf fLiio,L i Acu (a)(t/Li) 

Therefore by linearity of the Hilbert transform: 

ai(t/Li )) ~ L^ /2 g UJ (a)(t/L i ) 


□ 


Proposition B.8. The Hurst exponent H(f mp of g u (X(i)) satisfies the approx¬ 
imate relation: 


H a 'i r m 



if a < f3 + 3, a > 2 
otherwise 


(9) 


Proof. X(t) may be divided into a sum of long avalanches L S} i > L' which do 
not overlap and short avalanches. 


X(t) 


T & n . 


E 

L s i>L 


t — S 
Lad 



( 10 ) 
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Therefore by Proposition B.7 


ux(t))= £ /«(«., 

L s i>L' ' 


t — S 


E O 


i<U 


t — s 


( 11 ) 


The left hand term dominates whenever its variance is unbounded for large L c . 
This happens if r' < 3 which translates to a < (3 + 3 by Equation ([8]). Therefore 
the right hand term may be neglected and since the avalanches of the left hand 
term are separated in time, the envelope operator may be pulled under the sum 
so that: 


9u(X(t)) 


V rP /2 n 

/ j ^s,i 9^ 


i>L 


t — s 

Lei 


( 12 ) 


The same proof as for Proposition B.6 may then be applied to deriving H amp = 
^ — j. On the other hand if t' > 3 then the right hand term dominates, for 
large L' and has time-scale bounded by L' —thus in this case H amp = 1/2. □ 

We now study the DCCA correlation coefficients between amplitude en¬ 
velopes. 

Proposition B.9. Assuming that g u (a(t))dt converges, fora < /3 + 3 and 
a >2, p DC cA{n,g Ul {X{t)),g U 2 {X(t)) ->• 1 for n -)• oo. 

Proof. For large n, each DCCA window is longer than the largest avalanche. For 


a < (3 + 3 the variance in the left hand term of Equation (11) may be assumed 
to be non-overlapping since we assume a > 2 (Propostion B.5) we may neglect 


small avalanches and assume that a window contains only one avalanche at t' 
with length L. Then: 


J2^j(X){i) = 


0 


if t < t' 


i= 1 


Ll3/2 fo L ( a(t/L )) if t > t' + L 


But J 0 L g Ul ( a(t/L )) ~ g ( a(t/L )) up to a constant factor for large avalanches 
and assuming the integral converges. Thus the correlation tends to 1. □ 

Proposition B.10. Assuming that a(t)dt converges, for r < 3 and a > 2, 
PDCCA(n,g Ul (X(t)),g U2 (X(t)) -a 1 for n -> oo and ->• 0. 

Proof. For Wi,u )2 > L c then we have that g u .{L^a{ft/L)) ~ L 13 g Uj (L^a(t/L)). 
This is because, a(t/L) may be approximated by a delta function relative to 
these frequencies, which has white spectrum. Then the same proof techniques 
of Propositions B.8 and B.9 may be applied but with Equation replaced 
by: 


ux(t))= e O- 

L s ,i>L 



+ E L s,ifu 

L Sj i<L 
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For t > 3, the left hand term dominates, since it has unbounded variance in L c , 
which completes the proof. □ 

Proposition B.ll. For a < 2, E(p DC CA{n,g Ul (X(t)),g U 2 (X(t)))) -> 0 for 
n -> oo, Hf mp = 1/2 and 

Proof. The proof is identical to the proofs of Propositions |B.1| and |B.4| using 
the results of Proposition |B.6| □ 


B.l Summary of Results 

We term each region of parameters a meta-universality class which gives qual¬ 
itatively constant behaviour as quantified by Hurst exponents and DCCA cor¬ 
relation coefficients. There are four meta-universality classes MU1—MU4: 

MU1: a < 2: 


H r 


> 0.5 


Km V = 0-5 

Pdcca(ti, g Ul (X(t)),g U 2 (X(t))) -a 0 
MU2: a > 2 and a < (3 + 3 

H raw -' > 0*5 
tiam V > 0-5 


pDCCA{n, g Ul (X(t)),g U2 (X(t))) -> 1 


MU3: a > 2 and r < 3 


H raw ^ 0-3 

TTW _ q r 

amp yj ’ u 


PDCCA{n,g Ul (X(t)),g U 2 (X(t))) -> 1 as -a- 0 


MU4: a > 2 and r > 3 


Hr 


H: 


amp 


= 0.5 
= 0.5 


pDCCA{n, g Ul (X(t)),g U2 (X(t))) -a 0 


C Simulations 

In all simulations we model the average avalanche shape as a quadratic func¬ 
tion: b(t) = —4 (t — 1/2) 2 + 1. Moreover we set the variance swell identically 
so that b(t) = c{t). For the noise component we take the implementation of 
[261 . For the power-law cutoff sampling, we perform a density transformation 
of the uniform distribution. (In MATLAB x = rand(l,T) .*(1-L C )+L C ; x = 
l./(x.~ (1 ./(ck—1)))) 
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C.l Simulation 1 (from main text) 

We describe here how we obtain the results of Figure 2 of the main text. 

Lower panel: we filter white Noise with T = 40000 (using the method of BBl l 
to yield a process with spectrum scaling according to 1/w 19 . We then measure 
DCCA correlations between amplitudes in the frequency ranges [0.68,0.72] and 
[0.78,0.82] of half the sampling frequency, with n log-spaced between 20 and 
9000 and Hurst exponents of the same amplitudes using DFA with window sizes 
between 1000 and 9000. 

Top two panels: we generate two processes assuming CPLAD applied to 
Equation (1) of the main paper, generating avalanches as described above. We 
set q = 5, T = 40000 /? = 1 and a = 2.5,1.5 (top, middle resp.). We then apply 
the same analysis as in the lower panel. 

C.2 Simulation 2 

The aim of this simulation is to verify that long-range dependent Gaussian 
processes satisfy H^ mp = 0.5 and Pdcca = 0, where the time series are gener¬ 
ated as filtered Gaussian white noise processes. Thus this simulation tests our 
predictions for the PF model. 

To this end we simulate 100 long-range dependent Gaussian processes (H > 
0.5) (using the method of [26] ) of length 15000 time points. In each case the 
data are then filtered forward and backwards in two separate frequency bands 
(between 0.39 and 0.41 of the sampling frequency and 0.29 and 0.31 of the sam¬ 
pling frequency) with butterworth filters of order n. We then measure DCCA 
correlations between the Hilbert transforms of these signals and measure their 
Hurst exponents with DFA, in both cases using window lengths between 10 3 
and 10 4 . 

We repeat this setup for H = 0.9,1.4 and n = 2,4. 

The results are displayed in Figure [6] and show that, up to small sample 
effects, we expect H amp = 0.5 and zero cross correlations pdcca = 0 between 
frequency bands. 

C.3 Simulation 3 (from main text) 

In the first simulation, for each pair of exponents in the ranges a = 1.5,..., 6.5 
and /3 = 0.25,0.5,..., 3, we generate a sample path X(t) of length T = 300,000, 
with a cutoff at L c = 100,000, and number of superpositions q = 5. The 
first 100, 000 time points are discarded, to ensure stationarity. We then design 
Butterworth filters of order 2 between 0.29 and 0.31 and between 0.39 and 
0.41 of the sampling frequency. The data from X(t) are then filtered forwards 
and backwards (yielding effective filter order of 4) and the amplitude envelopes 
are calculated to yield Y\ (t) and Y 2 (t). We measure the Hurst exponent of 
£ ( X Ul ( t )), using DFA, and the DCCA correlation coefficients between £{X UJl ( t )) 
and £(X U2 (t)), setting n to log spaced values between 100 and 200000. This 
setup is repeated 100 times, and the results of the simulations are averaged. 
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Gaussian processes: H = 0.9 


Gaussian processes: H = 1.4 




Figure 6: The figure displays the results of the long-range analysis for a Gaussian 
process model. 

C.4 Simulation 4 

In the second simulation we check the prediction of H amp and H raw , setting set 
a = 2.5, q = 1, L c = 10 6 , with a burn in time of 10 5 time points (less time 
is required for convergence for this value of a), setting n to log spaced values 
between 7000 and 60000 for the estimation of H amp and between 100 and 7000 
for the estimation of H raw (according to where scaling regions were observed). 
The results are displayed in Figure [7] 

The quality of the H raw estimate is greater for small /3, whereas the quality 
of the H amp estimate is greater for larger (3. This discrepancy may be explained 
as follows: since X(t) has longer tails than £(X u (t)), the convergence of its 
moments is slow for large /3, thus the quality of the estimate decreases for large 
(3. On the other hand, the estimate of H raw requires a linear approximation 
to the non-linear transform given by the amplitude of the analytic signal: this 
approximation increases in quality for larger (3. 

C.5 Simulation 5 

The aim of this simulation is to verify the theory of Proposition |B.7| checking 
the heights of the filtered avalanches a u (t). We set (3 = 0.25, L = 2 10 ,... ,2 13 
and for each L considered, we simulate 10 4 avalanches of this length, and calcu¬ 
late the mean avalanche profile. We then log-regress the height of these profiles 
against log(L). The results are displayed in Figure [8] Close agreement is ob¬ 
served between the theoretical estimate (3' = (3/2 and the simulated results; the 
prediction improves for higher (3. 
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H amp ex P onents for a = 2.5; varying p 



Figure 7: Scaling of H amp and H raw (second simulation) compared to theory. 
The theoretical estimate is H amp ~ 2 — a/2 + /3' for (3' = 1. 


Scaling of filtered avalanche heights 



log(L) 


Figure 8: Comparison of the estimate /3/2 for the scaling of heights of filtered 
avalanches vs. simulation. 
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D Data Analysis 

Seven subjects participated in the study (1 female). The experimental pro¬ 
tocol was approved by the Institutional Review Board of the Charite Medical 
University, Berlin. EEG recordings were obtained at rest with subjects seated 
comfortably in a chair with their eyes open. Recordings were made of three ses¬ 
sions, each 5 minutes long so that each data set comprises roughly 15 minutes 
of data. EEG data were recorded with 96 Ag/AgCl electrodes, using BrainAmp 
amplifiers and BrainVision Recorder software (Brain Products GmbH, Munich, 
Germany). The signals were recorded in the 0.016-250 Hz frequency range at a 
1000Hz sampling frequency and subsequently subsampled to 200Hz. 

The data analytic steps taken on the EEG data were as follows. Outlier 
channels were rejected after visual inspection for abrupt shifts in voltage and 
poor signal quality. The data were then re-referenced according to the com¬ 
mon average. Spatial filters were computed on the data using Spatio-Spectral 
Decomposition (SSD) [27] . in order to extract components with pronounced al¬ 
pha oscillations. Spatial filters with poor signal quality or topographies were 
rejected. We then restricted our analysis to those filters displaying a peak in 
the alpha range; this step ensured a high signal quality with low levels of arti- 
factual activity. The fact that the spatial filters yield clear oscillatory signals 
ensured that the neuronal processes in the adjacent frequency ranges similarly 
originated from cortical areas relating to neuronal rather than artifactual activ¬ 
ity. For DFA and DCCA estimation we set n to log-spaced values between 1000 
and 25000. 


D.l Further Analyses 


In this section we repeat the analysis of real data reported in the main paper, but 
with a differing choice of spatial filter. Instead of taking SSD filters, which use 
the signal-to-noise ration in the alpha range to obtain filters corresponding to 
neuronal activity, we take Laplacian spatial filters (28j , which reduce redundancy 
between the data recorded at each electrode. All other parameters remain fixed. 
In this way, we show that the results we obtain do not depend on the alpha 
rhythm. 

The results are displayed in Figure[9]and confirm that the results obtained on 
real data using the SSD filters calibrated to the alpha rhythm may be obtained 
using filters which are independent of the alpha rhythm. 

To further check that results do not depend on the detrending degree d, we 
repeat this analysis with d = 2. The results are displayed in Figure [lO] and show 
qualitatively identical trends; most notable is that the results are less noisy 
for higher degree detrending, which we conjecture is due to a more thorough 
removal of trends for quadratic detrending. 

We also compute power-spectra of the first SSD component for each of the 
7 subjects, using Welch’s method (Hanning window of length 2 seconds, with 
0 overlap between windows and calculated at 256 points). The spectra are 


displayed in Figure 11 each spectrum clearly displays the presence of alpha/ 
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Figure 9: The figure displays results in an identical analysis as carried out for 
Figure 4 of the main paper, with the only difference being that we use Laplacian 
rather then SSD spatial filters. 

mu range oscillations at around 10Hz. 
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Figure 10: The figure displays results in an identical analysis as carried out for 
Figure 4 of the main paper, with the only difference being that we use Laplacian 
rather then SSD spatial filters. 
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PSD of 1st SSD: all subjects 
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Figure 11: The figure displays power-spectra from each of the 7 subjects con¬ 
sidered in the study. In case the power-spectrum is estimated on the first SSD 
component using Welch’s method. 

















